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Abstract 

We present the results of extensive computer simulations performed on 
solutions of monodisperse charged rod-like polyelectrolytes in the presence 
of trivalent counterions. To overcome energy barriers we used a combi- 
nation of parallel tempering and hybrid Monte Carlo techniques. Our 
results show that for small values of the electrostatic interaction the solu- 
tion mostly consists of dispersed single rods. The potential of mean force 
between the polyelectrolyte monomers yields an attractive interaction at 
short distances. For a range of larger values of the Bjerrum length, we find 
finite size polyelectrolyte bundles at thermodynamic equilibrium. Further 
increase of the Bjerrum length eventually leads to phase separation and 
precipitation. We discuss the origin of the observed thermodynamic sta- 
bility of the finite size aggregates. 

1 Introduction: 

Aggregation phenomena in charged macromolecular systems are well known 
and have been the focus of many investigations, see recent theoretical reviews 
12 El El ■ Here, we are interested in understanding the aggregation behav- 
ior of rod-like charged polyelectrolytes, such as DNA 0, F-actin|Hl d, and 
microtubules,^ that self-organize into bundles under the influence of only elec- 
trostatic interactions. Probably the most investigated systems are solutions of 
DNA in the presence of multivalent counterions ( see Refs. in jSl), which form 
a variety of different morphologies such as single or multiple toroids, or bundle- 
like aggregates. The origin of the attraction are short range ionic correlations of 
multivalent counterions, a phenomenon that cannot be explained on the level of 
mean-field theories like Poisson-Boltzmann ^|21E10- Simulations have already 
demonstrated this effect a long time ago |H1, but with the uprise of biological 
physics, it has regained interest in recent years. Many theoretical and compu- 
tational studies have dealt with the origin of the attraction, and we have an 
excellent qualitative understanding of it (sec Rcfs. in 01213101). 
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In a solution of many macromolecules, where multi-body effects are present, 
it has been difficult to obtain accurate enough analytical insight into the aggre- 
gate size distribution. The aggregate size must be determined by the competi- 
tion between the surface tension (due to short range attractions), the repulsive 
self-energy of the backbone charges, and the entropic degrees of freedom of the 
counterions. These effects are all coupled, and therefore pose a considerable chal- 
lenge to analytical approaches. For DNA, it has been argued |10[f51 [mil2| that 
the observed finite size of the aggregates is a kinetic effect based on the forma- 
tion of an energy barrier, and not a consequence of equilibrium thermodynamic 
properties of the system. Experiments performed with viral DNA support this 
conjecture at least for the investigated systems jj^l- However, recently Henle 
and Pincus argued that, depending on the actual parameters of the system, 
finite or infinite bundles should be possible in the presence of short-ranged at- 
tractions. Treating the system as consisting of sticky charged rods brings up 
the analogy to the splitting of a Rayleigh charged hydrophobic droplet, where 
the mean- field model has already been solved by Deserno ■ He observes that 
the droplet size is always finite, if the counterions cannot penetrate into the 
droplet, but that allowing counterion penetration leads either to finite or infi- 
nite droplet sizes, depending on the parameters. There have been relatively few 
simulations on bundle formations, notable exception being the simulations by 
Stevens ^1E1> showing the possibility that multivalent ions alone can lead to 
bundle formations, and the work of Borukhov et al. |18| . which focuses on two- 
rod systems bundled via short range mobile linker interactions. In the present 
work we investigate the bundle formation for a system of monodisperse charged 
semi-flexible polymers in the presence of trivalent counterions. This system was 
chosen since it is known to have short range attractions if the electrostatic in- 
teractions are sufficiently large. We investigate the aggregate size distribution 
as a function of interaction strength, and demonstrate clearly the existence of 
thermodynamically stable finite size aggregates. 

Note that there are also synthetic rod-like polyelectrolytes, such as poly(para- 
phenylene) oligomers jj^l i that form stable well defined cylindrical micelles due 
to hydrophobic interaction of their side chains (2UI 12 1| . Even though the dom- 
inant cause of attraction is of non-electrostatic origin, their stability diagram 
has some relation to the DNA-like systems under investigation which is clear if 
one considers an Ansatz like in Ref. |15j . 

In the first section we will present the model for our computational approach. 
In the following results and discussion section we will present the phase diagram 
for this system, and explain the origin of the stability of the finite size aggregates. 
We end with our conclusions. 

2 Model: 

In this study a coarse-grained representation of polyelectrolyte (PE) chains 
and trivalent counterions is used. All particles interact with purely repulsive 
Lennard- Jones interactions (Eqn. ^ with cut-off distance rcut/c = "2^1^ ^ in- 
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teraction strength clj = ksT, and the same diameter a. Each PE chain is 
composed of 30 negatively charged beads connected with FENE bonds (Eqn. [SJ 
with a spring constant of kp = TksT and cut-off distance Rp/cr ~ 2. The PE 
chains are semi-fiexible with a harmonic bending potential (Eqn. |3Jl of stiffness 
ke = IOO/cbT. We are interested in the dilute regime where the bundle-bundle 
contacts could be ignored. The PE monomer number density is fixed for all 
simulations as p/(J~^ = 7.5 10~^ . There are a total number of 61 PE chains 
and 10 trivalent counterions per chain in the cubic simulation box, where pe- 
riodic boundary conditions were used. The unscreened coulomb interactions 
(Eqn. 2J) are calculated using a version of the P3M method In Eqn. 21 
is the Bjerrum length and is used to control the strength of the electrostatic 
interactions. We used MD simulations with a time step of At — O.OOSr, where 
T is the usual Lennard-Jones time unit, using the Espresso paclcage|23|. 



ULj{r) ^ icLj ([^y^ - [^y) for r<re„t (1) 




Upir) = ^kpRl In 1^1 - ( — ) ) for r < Rp (2) 

Ue{r) = \kge^ (3) 

UciUj) = XsksT^ (4) 

The simulation of stiff-chain PEs with trivalent counterions possess a multi- 
tude of difficulties because of the long-range electrostatic interactions and energy 
barriers. In our previous study[2ni on bundles of hydrophobically modified PEs, 
we have shown that one can obtain a phase diagram by studying the stability of 
different size aggregates as a function of hydrophobic interaction strength and 
Bjerrum length. In that study, MD simulations with a Langevin thermostat 
were used to look at the stability of a preformed assembly of the PEs. 

Simulations of the current PE system without explicit hydrophobic inter- 
actions with the same approach as in Ref. [201 showed that, for low values of 
the Bjerrum length {Xb/c < 1-70), bundles are not stable. The initial bundle 
quickly falls apart. At equilibrium one occasionally observes short-lived dimers 
and trimers, but bigger aggregates are not observed. On the other hand for 
Ab/ct > 1.80 the situation changes drastically. If one starts from aggregates up 
to size A'' = 12 (iV is the number of PEs in the bundle), a fraction of PEs split 
up, but a finite size aggregate remains in solution. The size of the remaining 
aggregate depends on the size of the starting bundle. However, as the aggregate 
size is set to > 19 the bundles remain stable throughout the simulation. Be- 
yond Ab/(t = 2.0 independent of the bundle size, no PEs disintegrate from the 
start-up bundle. The reason for this interesting behavior could be explained by 
two different scenarios. By looking at this as a nucleation problem one can argue 
that for Ab/c = 1.8 — 2.0 the critical nucleus size is larger than 19 molecules. 
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Since the bundles we have studied are smaller than the critical nucleus size 
the bundles never fall apart completely or grow bigger, but one observes large 
fluctuations in the bundle size. We have tested this by studying an aggregate 
of size 25 with additional dispersed rods in solution for As/cr > 2.0. After an 
integration time of 5 10'*t no aggregate growth is observed, which leads us to 
conclude that this is not just a nucleation phenomenon. 

Another explanation for this behavior could be high energy barriers for the 
PEs to split up. This led us to believe that the Langevin thermostat is not 
efficient enough to overcome such energy barriers, which are highly dependent 
on Ab and aggregate size. The parallel tempering method (PT) |^ I^HI with 
Bjerrum length as the tempering parameter was chosen to overcome such energy 
barriers. We have set up 49 ensembles with Ab/(t ranging from 1.50 — 2.19. For 
the simulation of the individual ensembles, we have employed a hybrid Monte 
Carlo method {HMC)|2niE3- We integrated for 50 r for each HMC trial, which 
yields a good combination of acceptance ratio and efficient sampling of the con- 
figurational space. For each PT exchange attempt, we have done 9 HMC moves, 
which we call one cycle. The acceptance ratio of the HMC moves largely depends 
on the aggregate size distribution, and decreases rapidly for large bundle sizes. 
In order to test the efficiency of the PT+HMC method for the PE chains, we 
have conducted two separate simulations with 19 molecules each. The starting 
configurations were chosen as a single bundle for the first simulation and ran- 
domly distributed PEs for the second. After both of the PT+HMC simulations 
are equilibrated, an identical aggregate size distribution is observed. The re- 
sults presented in the following section are obtained from PT-I-HMC simulation 
with 61 PEs and neutralizing trivalent counterions. All ensembles are started 
from a configuration of randomly distributed rods. The set of ensembles have 
been equilibrated for up to 3500 cycles, and data is collected from an additional 
run of at least 2500 cycles, which took roughly three cpu years on Intel Xeon 
processors (2.4 GHz). 



3 Results and Discussion: 

We will start our analysis with the aggregate size distribution obtained from 
the PT-I-HMC simulation for the range of Bjerrum length values (Ab/ct = 
1.50 — 2.19). In order to define the aggregation state of the system, we as- 
sumed that any trivalent counterion closer than 3(t is condensed onto a PE. If 
two PEs share one or more condensed counterions they are assumed to be part 
of the same bundle. The number average {Nn — J2 ^ifi) and weighted average 
{Nw = fi/ ^Nifi) aggregate size (where fi is the fraction of aggregates 

of size Ni), as well as the polydispersity index (Nw/Nn) is shown in Fig. ^ 
In the weak electrostatic regime {Xb/ct 0) one observes only single PEs in 
solution, and therefore the polydispersity index is exactly one. For low Xb 
(up to Ab/c = 1.80) mostly single dispersed rods and occasionally short-lived 
dimers and trimers can be seen in solution, in agreement with our earlier ob- 
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Ab/(T, Bjerrum length 



Figure 1: The weighted {N^,) and number average {Nn) bundle size as a func- 
tion of Bjerrum length (As). For As/cr « 1.80 — 2.00 bundles of finite size are 
observed at thermal equilibrium. The Nw increases with As and the polydis- 
persity index has a maximum at As/cr w 1.95. 



servations. At such low values of As, even though the net charge of the rods is 
highly reduced by the condensed counterions, the repulsive interactions among 
PEs dominate the system behavior, and the polydispersity index is very close 
to one. Beyond As/f ~ 1-80 a gradual increase in the average aggregate size is 
observed. This clearly shows that finite size PE bundles exist at thermodynamic 
equilibrium for a range of Ab values. In this range the electrostatic interactions 
provide both the glue required for aggregation through short range correlations, 
and also the repulsive forces, which terminate aggregate growth. The entropic 
penalty still keeps counterions in solution and prevents complete phase separa- 
tion. Upon aggregation the polydispersity index rises above one, and shows a 
peak value of 3.7 around Ab/ct = 1.97. In the high As regime (As/cr > 2.1) the 
correlational attractions dominate, the entropic terms become negligible, and 
one expects the formation of an infinitely large bundle with a dilute solution 
of PEs. Therefore the polydispersity index should go to one again, which it 
does. Since we have only 61 PEs in our simulation box, once the average aggre- 
gate size reaches above 20, finite size effects dominate, and one cannot obtain a 
conclusive answer for the actual average aggregate size. 

One can observe the onset of aggregation by looking at the potential of mean 
force (PMF) among the centers of mass of the rods. However, the PMF of two 
rods still depends on the relative orientations ^2EH|- Instead, we have chosen 
to look at the PMF of the PE monomers (Fig. |2J), which is obtained by inverting 
their radial distribution function. We have performed an extra simulation with 
Ab/c = 1-0 to demonstrate the purely repulsive regime. Already for As/cr = 
1.5, the effective charge of the PEs is highly reduced and at short distances a 
weak attractive interaction (w ksT) is observed. As seen in Fig. ^this weak 
attraction is not sufficient to keep even two PEs together. The presence of 
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Figure 2: The potential of mean force among PE monomers. The purely 
repulsive interaction (Ab/ct — 1.0) is replaced with a weak attraction already 
for Xb/ct — 1.5. A small energy barrier is observed at a finite distance. As 
Xb is increased further the attraction becomes stronger and the energy barrier 
disappears, leading to the onset of aggregation into PE bundles. 



a weak energy barrier is also observed for A^/ct = 1.5. As one increases Xb 
further the energy barrier disappears, and the binding energy due to counterion 
correlations is significantly increased, leading to the formation of stable finite 
size aggregates. 

The distribution of PEs among different size aggregates at fixed Xb reveals 
that this is indeed a micellization phenomenon [521 (Fig. I^J. For Xb/ct — 1.84 
more than 60% of the PEs exist as dispersed single rods, whereas long-lived 
dimers and trimers exist in equilibrium with the single PEs. Upon further 
increase of (Ab/ct = 1.94) a peak for a finite aggregate size iV « 17 is observed. 
At Xb/ct — 1.99 the distribution broadens, and no single distinguishable peak 
remains. 

Next, we will look at the effective line charge density (/e) of different size 
bundles. For calculating /e, we assume that a bundle can be approximated as 
a rod- like charged object. The net charge of the bundle (PEs and condensed 
counterions) is divided by the length of the PEs to obtain /e. Note that these 
bundles have a finite diameter and small aspect ratio. Furthermore, the PEs 
are not perfectly aligned, but one observes large fluctuations along the bundle 
axis. In Fig. 21 we plot /e for aggregates of size 1, 2, 5, 10, and 23 as a function 
of As- For single PEs /e decreases hnearly up to As /a « 1.8. Beyond this value 
we observe a quick convergence of /e to a plateau value of w 0.1. If we look at 
/e obtained for a single PE simulated at the same density within the cell model, 
no such convergence is observed. The cell model simulations match with many- 
rod PT-I-HMC simulations only for low A^. This mismatch demonstrates the 
importance of many body effects, such as aggregation of PEs into bundles. For 
the many-rod simulation we see that beyond Xb/u > 1.8 most of the single PEs 
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Figure 3: The distribution of PEs among different size aggregates. Beyond the 
critical As a micelhzation type size distribution is observed, where monomers 
are in cquihbrium with finite size PE bundles. 
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Figure 4: Effective line charge density (/e) as a function of As for aggregates 
composed of N rods. At As ~ 1.8cr the average aggregate size strongly increases, 
which leads to a flattening of /g. For comparison we show fe of a single PE in 
the cell model, where multi-body effects are not present. In the insert the total 
fraction of condensed counterions is shown as a function of A^. 
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compensate a high fraction of their intrinsic charge via condensing counterions. 
As a result the repulsion among PEs is greatly diminished and the single PEs 
with the lowest line charge density aggregate to form bundles, leaving behind 
only the single PEs with higher fg. Upon aggregation of PEs, the trivalent 
counterions gain enthalpy, without diminishing the system entropy further more. 
Therefore /e for single PEs deviates from the cell model. For aggregates of 
size of 2, 5, 10, and 23 a similar behavior is observed. The reduction in /e 
continues until merging into even bigger aggregates becomes feasible, lowering 
the free energy even more. Note that large aggregate sizes are only observed 
beyond a certain Xb value. For example, for N = 5 the high fluctuations in /e 
around Ab/ct « 1.7 stem from the fact that these aggregates are very rare at 
such Xb values. The insert in Fig. 0] denotes the total fraction of condensed 
counterions over all PEs. This fraction increases linearly up to As/cr « 1.8. 
Beyond As/cr « 2.0 the charge compensation begins to saturate. Note that 
at this point more than 95% of the charge on the PEs are compensated, and 
only a small fraction of counterions are left in solution. Therefore the entropic 
penalty associated with condensing these remaining free counterions is extremely 
high. The deviation from the linear regime nicely coincides with the onset of 
aggregation. 

If we want to map our results to biomolecules such as DNA in a trivalent 
counterion solution we can do so by requiring that the strong coupling parameter 
S ^ for this system is the same (S w 80). For our model this is achieved for 
Ab ~ 1.3, where we do not yet observe aggregation. This is due to the short 
length of our rods. For short rods counterion condensation is much weaker |3U) 
due to end effects, and therefore one needs stronger electrostatic interactions to 
observe attraction and bundle formation. 

4 Conclusions: 

We investigated the aggregation behavior of semi-flexible polyelectrolytes with 
trivalent counterions. In dilute solutions the system shows a micellization type 
aggregate size distribution. At low values of the Bjerrum length single dispersed 
PEs are found, whereas at the other extreme of large couplings macroscopic 
phase separation is observed. Our most important finding is that finite size 
aggregates exist at thermodynamic equilibrium for a small window of Bjerrum 
length values. The aggregate size distribution of these PE bundles is rather 
broad. The finite size aggregates are a result of the correlation of the trivalent 
counterions at short distances, and repulsive electrostatic interactions at long 
distances. The entropic terms due to counterions, and end-effects due to the 
finite length of the PE chains also contribute to this delicate balance. The 
fraction of condensed counterions are found to be significantly lower than those 
for an infinite PE rod, which leads to an increased counterion entropy. 
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